#!/usr/bin/perl

use DBI;
use POSIX;
use strict;

# get our command line arguments
my $output_debug=0;
my $make_list=0;
my @rpfits_files;
my %reduction_parameters;
my %overrides;
$reduction_parameters{'robust_value'}=1;
$reduction_parameters{'self_calibration_loops'}=0;
for (my $i=0;$i<=$#ARGV;$i++){
    if ($ARGV[$i] eq "-D"){
	# level up debug printing
	$output_debug++;
    } elsif ($ARGV[$i] eq "-I"){
	# the name of the source to image
	$i++;
	$reduction_parameters{'imaging_source'}=$ARGV[$i];
    } elsif ($ARGV[$i] eq "-P"){
	# the name of the phase calibrator to use
	$i++;
	$reduction_parameters{'phase_calibrator'}=$ARGV[$i];
    } elsif ($ARGV[$i] eq "-F"){
	# the name of the flux calibrator to use
	$i++;
	$reduction_parameters{'flux_calibrator'}=$ARGV[$i];
    } elsif ($ARGV[$i] eq "-B"){
	# the name of the bandpass calibrator to use
	$i++;
	$reduction_parameters{'bandpass_calibrator'}=$ARGV[$i];
    } elsif ($ARGV[$i] eq "-L"){
	# the user wants a list of all the sources in this file
	$make_list=1;
	$output_debug++;
    } elsif ($ARGV[$i] eq "-f"){
	# only reduce a specific frequency
	$i++;
	$reduction_parameters{'reduce_frequency'}=$ARGV[$i];
    } elsif ($ARGV[$i]=~/^\-O(.*)/){
	# the user wants to override a restriction
	my $override=$1;
	if ($override eq "B"){
	    # will use whatever the user wants as the bandpass calibrator
	    $overrides{'bandpass_calibrator'}=1;
	} elsif ($override eq "D"){
	    # if we find that data is already existing, don't redo some
	    # data extraction steps
	    $overrides{'no_data_loading'}=1;
	} elsif ($override eq "R"){
	    # don't remove any data along the way
	    $overrides{'no_remove_data'}=1;
	} elsif ($override eq "FB"){
	    # don't do any automatic flagging of bandpass calibrator
	    $overrides{'no_autoflag_bandpass'}=1;
	} elsif ($override eq "FP"){
	    # don't do any automatic flagging of phase calibrator
	    $overrides{'no_autoflag_phase'}=1;
	} elsif ($override eq "FF"){
	    # don't do any automatic flagging of flux calibrator
	    $overrides{'no_autoflag_flux'}=1;
	} elsif ($override eq "C"){
	    # don't use mfclean, use regular clean
	    $overrides{'no_mfclean'}=1;
	}
    } elsif ($ARGV[$i] eq "-T"){
	# the user has specified a date
	$i++;
	$reduction_parameters{'date'}=$ARGV[$i];
    } elsif ($ARGV[$i] eq "-R"){
	# use robust imaging to minimise sidelobes
	$i++;
	$reduction_parameters{'robust_value'}=$ARGV[$i];
    } elsif ($ARGV[$i] eq "-S"){
	# use a self-calibration technique
	$i++;
	$reduction_parameters{'self_calibration_loops'}=$ARGV[$i];
    } else {
	# the name of an RPFITS file to use in the reduction
	push @rpfits_files,$ARGV[$i];
    }
}

# database schema
my %observation;
($observation{'source'},$observation{'ra'},$observation{'dec'},$observation{'n_if'},
 $observation{'if_freq'},$observation{'if_nchans'},$observation{'if_bw'},$observation{'time'},
 $observation{'arrayname'},$observation{'project_code'},$observation{'rpfits_file'},
 $observation{'cycle_time'},$observation{'offsource'},$observation{'onsource'},
 $observation{'scantype'},$observation{'calcode'},$observation{'stations'},
 $observation{'header_number'},$observation{'num_processed'},$observation{'processed_ids'},
 $observation{'id'},$observation{'az'},$observation{'el'},$observation{'attenmm'},
 $observation{'subrefl'},$observation{'corr_cfg'},$observation{'airtemp'},
 $observation{'chi'},$observation{'evector'},$observation{'jyperk'},
 $observation{'pntra'},$observation{'pntdec'},$observation{'relhumid'},
 $observation{'smonrms'},$observation{'wind'},$observation{'winddir'},
 $observation{'xtsys'},$observation{'ytsys'},$observation{'xyphase'},$observation{'xcaljy'},
 $observation{'ycaljy'},$observation{'xgtp'},$observation{'ygtp'},$observation{'xsdo'},
 $observation{'ysdo'},$observation{'pointcor'},$observation{'point_time'},
 $observation{'point_source'},$observation{'point_az'},$observation{'point_el'},
 $observation{'usable'})=
    (0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,
     31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50);

# database details
my $dsn="dbi:Pg:dbname=atca;host=mentok;port=5432";
my $user="postgres";
my $password="";

# Setup the database connection
my $dbh=DBI->connect($dsn,$user,$password,{ RaiseError => 0, AutoCommit => 1});

if ($output_debug>0){
    print "MESSAGE: imager starts\n";
}

# check we don't have conflicting arguments
if ($reduction_parameters{'date'} &&
    ($#rpfits_files>=0)){
    print "ERROR: date was specified along with a list of RPFITS files\n";
    exit;
}

# are we getting data from a specific date?
if ($reduction_parameters{'date'}){
    # get a list of all the RPFITS files from that date
    my $rpd_query="SELECT rpfits_file FROM observation WHERE time>='".
	$reduction_parameters{'date'}." 00:00:00' AND time<='".
	$reduction_parameters{'date'}." 23:59:59' ORDER BY time";
    if ($output_debug>1){
	print "MESSAGE: executing query: $rpd_query\n";
    }
    my $rpd=$dbh->prepare($rpd_query);
    $rpd->execute;
    while (my @rpd_response=$rpd->fetchrow){
	my $already_on=0;
	for (my $i=0;$i<=$#rpfits_files;$i++){
	    if ($rpfits_files[$i] eq $rpd_response[0]){
		$already_on=1;
		last;
	    }
	}
	if ($already_on==0){
	    push @rpfits_files,$rpd_response[0];
	    if ($output_debug>0){
		print "LIST: RPFITS file = ".$rpd_response[0]."\n";
	    }
	}
    }
    $rpd->finish;
}

my $nsource=1;
my @source_list;
my @all_sources;
my @all_ids;
my @all_starttimes;
my @all_endtimes;
my @all_ifs;
# which mode are we operating in?
if ($#rpfits_files>=0){
    # user wants us to make an image from a supplied set of files
    # make a list of all the sources
    my $lst_query="SELECT source,ra,dec,time,time+totaltime(maxcycles(onsource),".
	"maxcycles(offsource),cycle_time),if_freq,id FROM observation WHERE (";
    for (my $i=0;$i<=$#rpfits_files;$i++){
	if ($i!=0){
	    $lst_query.=" OR ";
	}
	$lst_query.="contains_text(rpfits_file,'".$rpfits_files[$i]."')=true";
    }
    $lst_query.=") ORDER BY time";
    if ($output_debug>1){
	print "MESSAGE: executing query: $lst_query\n";
    }
    my $lst=$dbh->prepare($lst_query);
    $lst->execute;
    while (my @lst_response=$lst->fetchrow){
	if ($make_list==1){
	    if ($output_debug>0){
		print "LIST: source $nsource (".$lst_response[3]." -> ".$lst_response[4]."): ".
		    $lst_response[0]." [".$lst_response[1]." ".$lst_response[2]."] @ ".
		    $lst_response[5]."\n";
	    }
	}
	@source_list=&add_source($lst_response[0],@source_list);
	push @all_sources,$lst_response[0];
	push @all_ids,$lst_response[6];
	push @all_starttimes,$lst_response[3];
	push @all_endtimes,$lst_response[4];
	push @all_ifs,$lst_response[5];
	$nsource++;
    }
    $lst->finish;
    
    if ($make_list==1){
	exit;
    }
}

# do some consistency checks
# have we been given a source to image?
if (!$reduction_parameters{'imaging_source'}){
    print "ERROR: no source was specified for imaging\n";
    exit;
}

# are the specified sources in the files?
for (my $i=0;$i<=$#source_list;$i++){
    if ($output_debug>0){
	print "MESSAGE: found source ".$source_list[$i]."\n";
    }
}
my $source_found=&source_present($reduction_parameters{'imaging_source'},
				 @source_list);
if ($source_found==0){
    print "ERROR: imaging source is not present in the RPFITS files\n";
    exit;
}
if ($reduction_parameters{'phase_calibrator'}){
    $source_found=&source_present($reduction_parameters{'phase_calibrator'},
				  @source_list);
    if ($source_found==0){
	print "ERROR: phase calibrator source is not present in the RPFITS files\n";
	exit;
    }
}
if ($reduction_parameters{'flux_calibrator'}){
    $source_found=&source_present($reduction_parameters{'flux_calibrator'},
				  @source_list);
    if ($source_found==0){
	print "ERROR: flux calibrator source is not present in the RPFITS files\n";
	exit;
    }
}
if ($reduction_parameters{'bandpass_calibrator'}){
    $source_found=&source_present($reduction_parameters{'bandpass_calibrator'},
				  @source_list);
    if ($source_found==0){
	print "ERROR: bandpass calibrator source is not present in the RPFITS files\n";
	exit;
    }
}    

# check that the specified phase, flux and bandpass calibrators are phase, flux
# and bandpass calibrators
if ($reduction_parameters{'phase_calibrator'}){
    my $pc_type=&calibrator_type($dbh,$reduction_parameters{'phase_calibrator'});
    if ($pc_type eq ""){
	print "ERROR: phase calibrator is not a known calibrator\n";
	exit;
    }
}
if ($reduction_parameters{'flux_calibrator'}){
    my $pc_type=&calibrator_type($dbh,$reduction_parameters{'flux_calibrator'});
    if ($pc_type ne "f"){
	print "ERROR: flux calibrator is not a known flux calibrator\n";
	exit;
    }
}
if ($reduction_parameters{'bandpass_calibrator'}){
    my $pc_type=&calibrator_type($dbh,$reduction_parameters{'bandpass_calibrator'});
    if ($pc_type eq ""){
	if (!$overrides{'bandpass_calibrator'}){
	    print "ERROR: bandpass calibrator is not a known calibrator\n";
	    exit;
	} else {
	    if ($output_debug>0){
		print "WARNING: bandpass calibrator ".$reduction_parameters{'bandpass_calibrator'}.
		    " is likely unsuitable; proceeding due to user override\n";
	    }
	}
    }
}

if ($output_debug>0){
    print "MESSAGE: imaging source ".$reduction_parameters{'imaging_source'}."\n";
}

# what frequencies are we imaging?
my @unique_imaging_ifs;
for (my $i=0;$i<=$#all_sources;$i++){
    if ($all_sources[$i] eq $reduction_parameters{'imaging_source'}){
	my $already_printed=0;
	for (my $j=0;$j<=$#unique_imaging_ifs;$j++){
	    if ($unique_imaging_ifs[$j] eq $all_ifs[$i]){
		$already_printed=1;
		last;
	    }
	}
	if ($already_printed==0){
	    push @unique_imaging_ifs,$all_ifs[$i];
	    print "MESSAGE: imaging source found with frequency configuration: ".
		$all_ifs[$i]."\n";
	}
    }
}
if (!$reduction_parameters{'reduce_frequency'}){
    if ($#unique_imaging_ifs>0){
	# too many IF configurations, what would we do?!
	print "ERROR: too many IF configurations; you must specify which to image\n";
	exit;
    }
    # choose the frequency of the first IF as the reduction frequency
    my @if_freqs=&postgresql_array_get($unique_imaging_ifs[0]);
    $reduction_parameters{'reduce_frequency'}=$if_freqs[0];
    $reduction_parameters{'reduce_if'}=1;
} else {
    # check that the frequency supplied is one of the available frequencies
    my $found_frequency=0;
    for (my $i=0;$i<=$#unique_imaging_ifs;$i++){
	my @if_freqs=&postgresql_array_get($unique_imaging_ifs[0]);
	for (my $j=0;$j<=$#if_freqs;$j++){
	    if ($reduction_parameters{'reduce_frequency'}==$if_freqs[$j]){
		$reduction_parameters{'reduce_if'}=$j+1;
		$found_frequency=1;
		last;
	    }
	}
	if ($found_frequency==1){
	    last;
	}
    }
    if ($found_frequency==0){
	# didn't find that frequency
	print "ERROR: could not identify any observations with frequency ".
	    $reduction_parameters{'reduce_frequency'}."\n";
	exit;
    }
}
if ($output_debug>0){
    print "MESSAGE: will image frequency ".$reduction_parameters{'reduce_frequency'}.
	" MHz (IF #".$reduction_parameters{'reduce_if'}.")\n";
}
# is the bandpass calibrator appropriate for this frequency
if ($reduction_parameters{'bandpass_calibrator'}){
    my $isbandpass=&is_bandpass_calibrator($dbh,$reduction_parameters{'bandpass_calibrator'},
					   $reduction_parameters{'reduce_frequency'});
    if ($isbandpass==0){
	if (!$overrides{'bandpass_calibrator'}){
	    # no good
	    print "ERROR: cannot use ".$reduction_parameters{'bandpass_calibrator'}." as a ".
		"bandpass calibrator at ".$reduction_parameters{'reduce_frequency'}." MHz\n";
	    exit;
	} else {
	    if ($output_debug>0){
		print "WARNING: bandpass calibrator ".$reduction_parameters{'bandpass_calibrator'}.
		    " is likely unsuitable; proceeding due to user override\n";
	    }
	}
    } else {
	if ($output_debug>0){
	    print "MESSAGE: using ".$reduction_parameters{'bandpass_calibrator'}." as the ".
		"bandpass calibrator\n";
	}
    }
}

# figure out the best sources to use if the user didn't specify them
# the phase calibrator
if (!$reduction_parameters{'phase_calibrator'}){
    if ($output_debug>0){
	print "MESSAGE: attempting to determine most suitable phase calibrator\n";
    }
    my @before_calibrators;
    my @after_calibrators;
    for (my $i=0;$i<=$#all_sources;$i++){
	if ($all_sources[$i] eq $reduction_parameters{'imaging_source'}){
	    my ($cal_before,$cal_after)=&nearest_source($dbh,$all_starttimes[$i],
							$all_endtimes[$i],$all_sources[$i]);
	    push @before_calibrators,$cal_before;
	    push @after_calibrators,$cal_after;
	}
    }
    my @count_source;
    my @source_count;
    for (my $i=0;$i<=$#before_calibrators;$i++){
	my $new_source=1;
	for (my $j=0;$j<=$#count_source;$j++){
	    if ($count_source[$j] eq $before_calibrators[$i]){
		$new_source=0;
		$source_count[$j]++;
		last;
	    }
	}
	if ($new_source==1){
	    push @count_source,$before_calibrators[$i];
	    push @source_count,1;
	}
    }
    for (my $i=0;$i<=$#after_calibrators;$i++){
	my $new_source=1;
	for (my $j=0;$j<=$#count_source;$j++){
	    if ($count_source[$j] eq $after_calibrators[$i]){
		$new_source=0;
		$source_count[$j]++;
		last;
	    }
	}
	if ($new_source==1){
	    push @count_source,$after_calibrators[$i];
	    push @source_count,1;
	}
    }
    my $max_count=0;
    my $max_position=-1;
    
    for (my $i=0;$i<=$#count_source;$i++){
	if ($source_count[$i]>$max_count){
	    $max_count=$source_count[$i];
	    $max_position=$i;
	}
    }
    $reduction_parameters{'phase_calibrator'}=$count_source[$max_position];
    if ($output_debug>0){
	print "MESSAGE: recommend using ".$reduction_parameters{'phase_calibrator'}.
	    " as the phase calibrator\n";
    }
} else {
    if ($output_debug>0){
	print "MESSAGE: using ".$reduction_parameters{'phase_calibrator'}." as the ".
	    "phase calibrator\n";
    }
}

# the flux calibrator
if (!$reduction_parameters{'flux_calibrator'}){
    if ($output_debug>0){
	print "MESSAGE: attempting to determine most suitable flux calibrator\n";
    }
    my $first_source=-1;
    for (my $i=0;$i<=$#all_sources;$i++){
	if ($all_sources[$i] eq $reduction_parameters{'imaging_source'}){
	    $first_source=$i;
	    last;
	}
    }
    my @source_ifs=&postgresql_array_get($all_ifs[$first_source]);
    my @near_events=&find_nearest_events($dbh,$all_starttimes[$first_source],
					 @source_ifs);
    my @fluxcal_ids=&find_nearest_fluxcals($dbh,$all_starttimes[$first_source],
					   @near_events,@source_ifs);
    if ($#fluxcal_ids>=0){
	my @id_response=&get_id_info($dbh,$fluxcal_ids[0]);
	$reduction_parameters{'flux_calibrator'}=$id_response[$observation{'name'}];
	if ($output_debug>0){
	    print "MESSAGE: recommend using ".$reduction_parameters{'flux_calibrator'}.
		" as the flux calibrator\n";
	}
    } else {
	if ($output_debug>0){
	    print "ERROR: unable to find suitable flux calibrator!\n";
	}
	exit;
    }
    # check that the flux calibrator is part of the observations
    $source_found=&source_present($reduction_parameters{'flux_calibrator'},
				  @source_list);
    if ($source_found==0){
	if ($output_debug>0){
	    print "MESSAGE: flux calibrator cannot be obtained from supplied files\n";
	}
	my @required_files;
	for (my $i=0;$i<=$#fluxcal_ids;$i++){
	    my @id_response=&get_id_info($dbh,$fluxcal_ids[$i]);
	    my $already_required=0;
	    for (my $j=0;$j<=$#required_files;$j++){
		if ($required_files[$j] eq $id_response[$observation{'rpfits_file'}]){
		    $already_required=1;
		    last;
		}
	    }
	    if ($already_required==0){
		push @required_files,$id_response[$observation{'rpfits_file'}];
	    }
	}
	for (my $i=0;$i<=$#required_files;$i++){
	    if ($output_debug>0){
		my $this_file=$required_files[$i];
		$this_file=~s/^(.*)\/([^\/]*)$/$2/;
		print "MESSAGE: require file ".$this_file."\n";
		push @rpfits_files,$required_files[$i];
	    }
	}
#	if ($output_debug>0){
#	    print "ERROR: cannot continue without flux calibrator\n";
#	}
#	exit;
    }
} else {
    if ($output_debug>0){
	print "MESSAGE: using ".$reduction_parameters{'flux_calibrator'}." as the ".
	    "flux calibrator\n";
    }
}

# the bandpass calibrator
if (!$reduction_parameters{'bandpass_calibrator'}){
    if ($output_debug>0){
	print "MESSAGE: attempting to determine most suitable bandpass calibrator\n";
    }
    for (my $i=0;$i<=$#source_list;$i++){
	my $isbandpass=&is_bandpass_calibrator($dbh,$source_list[$i],
					       $reduction_parameters{'reduce_frequency'});
	if ($isbandpass==1){
	    # we've found a suitable bandpass calibrator
	    $reduction_parameters{'bandpass_calibrator'}=$source_list[$i];
	    last;
	}
    }
    if ($reduction_parameters{'bandpass_calibrator'}){
	if ($output_debug>0){
	    print "MESSAGE: using ".$reduction_parameters{'bandpass_calibrator'}." as the".
		" bandpass calibrator\n";
	}
    } else {
	print "ERROR: unable to find a suitable bandpass calibrator!\n";
	exit;
    }
}

# we now attempt to find the RPFITS files
my @rpfits_required;
for (my $i=0;$i<=$#rpfits_files;$i++){
    $rpfits_files[$i]=~s/^.*\/([^\/]*)$/$1/;
    if (!-e $rpfits_files[$i]){
	push @rpfits_required,splice @rpfits_files,$i,1;
	$i--;
    }
}
# check which of the required files actually has the sources we
# require
for (my $i=0;$i<=$#rpfits_required;$i++){
    my $file_required=0;
    my $flr_query="SELECT source FROM observation WHERE contains_text(rpfits_file,'".
	$rpfits_required[$i]."')=TRUE AND if_freq='".$unique_imaging_ifs[0]."'";
    if ($output_debug>1){
	print "QUERY: $flr_query\n";
    }
    my $flr=$dbh->prepare($flr_query);
    $flr->execute;
    while (my @flr_response=$flr->fetchrow){
	if (($flr_response[0] eq $reduction_parameters{'phase_calibrator'})||
	    ($flr_response[0] eq $reduction_parameters{'flux_calibrator'})||
	    ($flr_response[0] eq $reduction_parameters{'bandpass_calibrator'})||
	    ($flr_response[0] eq $reduction_parameters{'imaging_source'})){
	    # this file is required
	    $file_required=1;
	    last;
	}
    }
    $flr->finish;
    if ($file_required==0){
	if ($output_debug>0){
	    print "MESSAGE: RPFITS file ".$rpfits_required[$i]." is not required for reduction\n";
	}
	splice @rpfits_required,$i,1;
	$i--;
    }
}
for (my $i=0;$i<=$#rpfits_required;$i++){
    if ($output_debug>0){
	print "MESSAGE: attempting to find RPFITS file ".$rpfits_required[$i]."\n";
    }
    $rpfits_required[$i]=&obtain_rpfits($rpfits_required[$i],".");
    if (-e $rpfits_required[$i]){
	# the file was copied successfully
	push @rpfits_files,splice @rpfits_required,$i,1;
	$i--;
    }
}
# did we get all the files?
if ($#rpfits_required>=0){
    print "ERROR: unable to locate all required RPFITS files\n";
    exit;
}

# we have our datasets, we know which calibrators to use, so let's get to it


# MIRIAD will output all files in lowercase, so we convert all our source names
# to lowercase
$reduction_parameters{'bandpass_calibrator'}=~tr/A-Z/a-z/;
$reduction_parameters{'flux_calibrator'}=~tr/A-Z/a-z/;
$reduction_parameters{'phase_calibrator'}=~tr/A-Z/a-z/;
$reduction_parameters{'imaging_source'}=~tr/A-Z/a-z/;

# first, we load the data into MIRIAD
if ($output_debug>0){
    print "MESSAGE: loading data into MIRIAD format\n";
}
my $uvdata=&miriad_load_data($reduction_parameters{'reduce_frequency'},
			     $reduction_parameters{'reduce_if'},
			     $overrides{'no_data_loading'},@rpfits_files);

# then we flag based on calibration events that we know about
if ($output_debug>0){
    print "MESSAGE: flagging setup data\n";
}
&event_flag($dbh,$uvdata,@rpfits_files);
		  
# then we split the data
if ($output_debug>0){
    print "MESSAGE: splitting dataset\n";
}
&miriad_split_data($uvdata,$overrides{'no_remove_data'});

# determine the reference antenna
$reduction_parameters{'reference_antenna'}=
    &determine_reference_antenna($dbh,@rpfits_files);
if ($output_debug>0){
    print "MESSAGE: using antenna ".$reduction_parameters{'reference_antenna'}.
	" as the reference\n";
}

# first step in reduction is to calibrate the bandpass calibrator
if ($output_debug>0){
    print "MESSAGE: determining bandpass solution\n";
}
if ($output_debug>0){
    print "MESSAGE: calibrating...\n";
}
&miriad_calibrate($reduction_parameters{'bandpass_calibrator'},
		  $reduction_parameters{'reduce_frequency'},
		  $reduction_parameters{'reference_antenna'},
		  1,0,0);

# and now we automatically flag bad data from this calibrator
if (!$overrides{'no_autoflag_bandpass'}){
    if ($output_debug>0){
	print "MESSAGE: flagging...\n";
    }
    &autoflag($reduction_parameters{'bandpass_calibrator'},
	      $reduction_parameters{'reduce_frequency'});


    # redo the calibration now with the flagged dataset
    if ($output_debug>0){
	print "MESSAGE: calibrating...\n";
    }
    &miriad_calibrate($reduction_parameters{'bandpass_calibrator'},
		      $reduction_parameters{'reduce_frequency'},
		      $reduction_parameters{'reference_antenna'},
		      1,0,0);
}

# next we copy the bandpass solution to the flux and phase calibrators
if ($output_debug>0){
    print "MESSAGE: copying bandpass solution to flux and phase calibrators\n";
}
if ($reduction_parameters{'bandpass_calibrator'} ne
    $reduction_parameters{'flux_calibrator'}){
    &miriad_copy_calibration($reduction_parameters{'bandpass_calibrator'},
			     $reduction_parameters{'flux_calibrator'},
			     $reduction_parameters{'reduce_frequency'});
}

if ($reduction_parameters{'bandpass_calibrator'} ne
    $reduction_parameters{'phase_calibrator'}){
    &miriad_copy_calibration($reduction_parameters{'bandpass_calibrator'},
			     $reduction_parameters{'phase_calibrator'},
			     $reduction_parameters{'reduce_frequency'});
}

# and calibrate both the flux and phase calibrators
if ($output_debug>0){
    print "MESSAGE: calibrating flux and phase calibrators\n";
}
&miriad_calibrate($reduction_parameters{'flux_calibrator'},
		  $reduction_parameters{'reduce_frequency'},
		  $reduction_parameters{'reference_antenna'},
		  0,1,1);
&miriad_calibrate($reduction_parameters{'phase_calibrator'},
		  $reduction_parameters{'reduce_frequency'},
		  $reduction_parameters{'reference_antenna'},
		  0,1,0);

# now flag these calibrators
if (!$overrides{'no_autoflag_flux'}){
    if ($output_debug>0){
	print "MESSAGE: automatically flagging flux calibrator\n";
    }
    &autoflag($reduction_parameters{'flux_calibrator'},
	      $reduction_parameters{'reduce_frequency'});
}
if (!$overrides{'no_autoflag_phase'}){
    if ($output_debug>0){
	print "MESSAGE: automatically flagging phase calibrator\n";
    }
    &autoflag($reduction_parameters{'phase_calibrator'},
	      $reduction_parameters{'reduce_frequency'});
}

# and recalibrate
if (!$overrides{'no_autoflag_flux'}){
    if ($output_debug>0){
	print "MESSAGE: calibrating flux calibrator\n";
    }
    &miriad_calibrate($reduction_parameters{'flux_calibrator'},
		      $reduction_parameters{'reduce_frequency'},
		      $reduction_parameters{'reference_antenna'},
		      0,1,1);
}
if (!$overrides{'no_autoflag_phase'}){
    if ($output_debug>0){
	print "MESSAGE: calibrating phase calibrator\n";
    }
    &miriad_calibrate($reduction_parameters{'phase_calibrator'},
		      $reduction_parameters{'reduce_frequency'},
		      $reduction_parameters{'reference_antenna'},
		      0,1,0);
}

# now we set the correct flux scaling
if ($output_debug>0){
    print "MESSAGE: correcting flux scale and bandpasses for phase calibrator\n";
}
&miriad_fluxscale($reduction_parameters{'flux_calibrator'},
		  $reduction_parameters{'phase_calibrator'},
		  $reduction_parameters{'reduce_frequency'});

# copy the final calibration to the imaging source
if ($output_debug>0){
    print "MESSAGE: copying gain and bandpass calibration to imaging source\n";
}
if ($reduction_parameters{'phase_calibrator'} ne
    $reduction_parameters{'imaging_source'}){
    &miriad_copy_calibration($reduction_parameters{'phase_calibrator'},
			     $reduction_parameters{'imaging_source'},
			     $reduction_parameters{'reduce_frequency'});
}

# and apply it
if ($output_debug>0){
    print "MESSAGE: applying imaging source calibration solution\n";
}
&miriad_apply_calibration($reduction_parameters{'imaging_source'},
			  $reduction_parameters{'reduce_frequency'});

# calibration complete, so it's imaging time
my ($ra_axis,$dec_axis)=(-1,-1);
my $imaging_completed=0;
my $image_image;
while ($imaging_completed==0){
    if ($output_debug>0){
	print "MESSAGE: creating dirty map\n";
    }
    my ($image_map,$image_beam,$image_rms,$image_xpixels,$image_ypixels,$imagetripled)=
	&miriad_make_image($reduction_parameters{'imaging_source'},$ra_axis,$dec_axis,
			   $reduction_parameters{'reduce_frequency'},
			   $reduction_parameters{'robust_value'});
    
    # now deconvolve
    if ($output_debug>0){
	print "MESSAGE: deconvolving image\n";
    }
    if ($overrides{'no_mfclean'}==1){
	$image_xpixels=-1;
	$image_ypixels=-1;
    }
    my $image_model=&miriad_deconvolve_image($reduction_parameters{'imaging_source'},
					     $reduction_parameters{'reduce_frequency'},
					     $image_map,$image_beam,$image_rms,
					     $image_xpixels,$image_ypixels);
    
    # and restore the image
    if ($output_debug>0){
	print "MESSAGE: restoring image\n";
    }
    $image_image;
    ($image_image,$ra_axis,$dec_axis)=
	&miriad_restore_image($reduction_parameters{'imaging_source'},
			      $reduction_parameters{'reduce_frequency'},
			      $image_map,$image_beam,$image_model);

    # now decide whether the image has been made correctly
    if ((($image_xpixels!=-1)&&($image_ypixels!=-1))||
	($overrides{'no_mfclean'}==1)){
	$imaging_completed=1;
    }

    # make a measurement of the image size
    my $max_flux;
    ($ra_axis,$dec_axis,$max_flux)=&miriad_image_size($image_image);
    if ($imagetripled==1){
	$ra_axis/=3;
	$dec_axis/=3;
    }

    # are we self-calibrating?
    if ($reduction_parameters{'self_calibration_loops'}>0){
	$imaging_completed=0;
	&miriad_self_calibrate($image_model,$reduction_parameters{'imaging_source'},
			       $reduction_parameters{'reduce_frequency'},$max_flux,
			       $reduction_parameters{'reference_antenna'});
	$reduction_parameters{'self_calibration_loops'}--;
    }	
}

# we're out, so we've finished making the image!
if ($output_debug>0){
    print "MESSAGE: imaging completed, image is in ".$image_image."\n";
}

$dbh->disconnect;

exit;

sub postgresql_array_get {
    my ($arraystring)=@_;
    # turns a postgresql array string into a perl array
    if ($output_debug>2){
	print "attempting to turn $arraystring into array\n";
    }
    $arraystring=~s/[\{\}]//g;
    my @els=split(/\,/,$arraystring);
    return @els;
}

sub find_nearest_events {
    my ($dbh,$time,@ifs)=@_;

    my $time_difference="event_time-'".$time."'";
    my ($nearest_before,$nearest_after);
    my $evf_query="SELECT * FROM events WHERE (".
	"(event_type='acal' AND if_freq='".
	&array_to_text(@ifs)."')) AND EXTRACT(epoch FROM (".$time_difference."))<0".
	" ORDER BY abs(EXTRACT(epoch FROM (".$time_difference."))) LIMIT 1";
    if ($output_debug>1){
	print "executing query: $evf_query\n";
    }
    my $evf=$dbh->prepare($evf_query);
    $evf->execute;

    while(my @evf_response=$evf->fetchrow){
	$nearest_before=$evf_response[1];
    }
    $evf->finish;

    my $evf_query="SELECT * FROM events WHERE (".
	"(event_type='acal' AND if_freq='".
	&array_to_text(@ifs)."')) AND EXTRACT(epoch FROM (".$time_difference."))>0".
	" ORDER BY abs(EXTRACT(epoch FROM (".$time_difference."))) LIMIT 1";
    if ($output_debug>1){
	print "executing query: $evf_query\n";
    }
    my $evf=$dbh->prepare($evf_query);
    $evf->execute;

    while(my @evf_response=$evf->fetchrow){
	$nearest_after=$evf_response[1];
    }
    $evf->finish;

    return ($nearest_before,$nearest_after);
}

sub miriad_load_data {
    my ($frequency,$if_number,$no_reload_data,@rpfits_files)=@_;

    my $uv_dataset="loaded_data.uv";
    # remove any current dataset from the directory
    if ((-e $uv_dataset)&&($no_reload_data!=1)){
	my $remove_command="rm -rf ".$uv_dataset;
	if ($output_debug>1){
	    print "MESSAGE: executing ".$remove_command."\n";
	}
	system $remove_command;
    }
	
    if (($no_reload_data!=1)||(!-e $uv_dataset)){
	my $atlod_command="atlod in=";
	for (my $i=0;$i<=$#rpfits_files;$i++){
	    if ($i>0){
		$atlod_command.=",";
	    }
	    $atlod_command.=$rpfits_files[$i];
	}
	$atlod_command.=" out=".$uv_dataset;
	$atlod_command.=" options=birdie,rfiflag,noauto,xycorr,opcorr ";
	$atlod_command.="ifsel=".$if_number;
	
	if ($output_debug>1){
	    print "MESSAGE: executing ".$atlod_command."\n";
	}
	my @loaddata_output=&run_command($atlod_command);
    
	&print_output_lines("MIRIAD",@loaddata_output);
    }

    return $uv_dataset;
}

sub run_command {
    my ($command)=@_;

    my @command_output;
    open(IN,"-|")||exec $command." 2>&1";
    while(<IN>){
	chomp;
	push @command_output,$_;
    }
    close(IN);

    return @command_output;
}

sub find_nearest_fluxcals {
    my ($dbh,$time,$nearevent_before,$nearevent_after,@ifs)=@_;

    my $unixtime_reference=&time_to_unixtime($time);
    my @fluxcal_ids;
    my @fluxcal_sources;
    my @fluxcal_maxonsources;
    
    # find an observation with the correct IF configuration nearest to our
    # observation time that also falls between the known nearby event times
    my $endtime="time+totaltime(maxcycles(onsource),maxcycles(offsource),cycle_time)";
    my $fnd_query="SELECT observation.id,observation.source,observation.onsource FROM ".
	"observation,calibrators WHERE observation.source=calibrators.name AND ".
	"calibrators.type='f'";
    if ($nearevent_before ne ""){
	$fnd_query.=" AND (".$endtime.")>'".$nearevent_before."'";
    }
    if ($nearevent_after ne ""){
	$fnd_query.=" AND (".$endtime.")<'".$nearevent_after."'";
    }
    $fnd_query.=" AND if_freq='".&array_to_text(@ifs)."'".
	" ORDER BY abs(EXTRACT(epoch FROM ((".$endtime.")-'".$time."')))";
    if ($output_debug>1){
	print "executing query: $fnd_query\n";
    }
    my $fnd=$dbh->prepare($fnd_query);
    $fnd->execute;
    while (my @fnd_response=$fnd->fetchrow){
	push @fluxcal_ids,$fnd_response[0];
	push @fluxcal_sources,$fnd_response[1];
	my @onsource_values=&postgresql_array_get($fnd_response[2]);
	my ($min_ons,$max_ons)=&minmax(@onsource_values);
	push @fluxcal_maxonsources,$max_ons;
    }
    $fnd->finish;

    # figure out which flux calibrator has the most time onsource
    my @fluxcal_found_source;
    my @fluxcal_found_onsource;
    for (my $i=0;$i<=$#fluxcal_sources;$i++){
	my $already_found=-1;
	for (my $j=0;$j<=$#fluxcal_found_source;$j++){
	    if ($fluxcal_found_source[$j] eq $fluxcal_sources[$i]){
		$already_found=$j;
		last;
	    }
	}
	if ($already_found==-1){
	    push @fluxcal_found_source,$fluxcal_sources[$i];
	    push @fluxcal_found_onsource,$fluxcal_maxonsources[$i];
	} else {
	    $fluxcal_found_onsource[$already_found]+=$fluxcal_maxonsources[$i];
	}
    }
    my ($best_source,$best_onsource);
    $best_source=$fluxcal_found_source[0];
    $best_onsource=$fluxcal_found_onsource[0];
    for (my $i=1;$i<=$#fluxcal_found_source;$i++){
	if ($fluxcal_found_onsource[$i]>$best_onsource){
	    $best_onsource=$fluxcal_found_onsource[$i];
	    $best_source=$fluxcal_found_source[$i];
	}
    }

    # now remove all ids that don't correspond to the best source
    for (my $i=0;$i<=$#fluxcal_ids;$i++){
	if ($fluxcal_sources[$i] ne $best_source){
	    splice @fluxcal_ids,$i,1;
	    splice @fluxcal_sources,$i,1;
	    splice @fluxcal_maxonsources,$i,1;
	    $i--;
	}
    }

    return (@fluxcal_ids);
}

sub minmax {
    # get the minimum and maximum values in an array
    my @array=@_;

    my $minval=$array[0];
    my $maxval=$array[0];
    for (my $i=0;$i<=$#array;$i++){
	if ($array[$i]<$minval){
	    $minval=$array[$i];
	}
	if ($array[$i]>$maxval){
	    $maxval=$array[$i];
	}
    }
    return ($minval,$maxval);
}

sub nearest_source {
    my ($dbh,$start_time,$end_time,$source)=@_;
    
    my ($source_before,$source_after)=("","");

    my $nrp_query="SELECT source FROM observation,calibrators WHERE EXTRACT(epoch FROM (".
	"time-'".$start_time."'))<0 AND source!='".$source."' AND name=source ".
	"ORDER BY EXTRACT(epoch FROM (time-'".$start_time."')) DESC LIMIT 1";
    if ($output_debug>1){
	print "executing query: $nrp_query\n";
    }
    my $nrp=$dbh->prepare($nrp_query);
    $nrp->execute;
    while (my @nrp_response=$nrp->fetchrow){
	$source_before=$nrp_response[0];
	last;
    }
    $nrp->finish;

    my $nrl_query="SELECT source FROM observation,calibrators WHERE EXTRACT(epoch FROM (".
	"time-'".$end_time."'))>0 AND source!='".$source."' AND name=source ".
	"ORDER BY EXTRACT(epoch FROM (time-'".$end_time."')) LIMIT 1";
    if ($output_debug>1){
	print "executing query: $nrl_query\n";
    }
    my $nrl=$dbh->prepare($nrl_query);
    $nrl->execute;
    while (my @nrl_response=$nrl->fetchrow){
	$source_after=$nrl_response[0];
	last;
    }
    $nrl->finish;

    return ($source_before,$source_after);
}

sub calibrator_type {
    my ($dbh,$source)=@_;
    
    my $type="";
    my $chk_query="SELECT type FROM calibrators WHERE name='".$source."'";
    if ($output_debug>1){
	print "executing query: $chk_query\n";
    }
    my $chk=$dbh->prepare($chk_query);
    $chk->execute;
    while (my @chk_response=$chk->fetchrow){
	$type=$chk_response[0];
    }
    $chk->finish;

    return $type;
}

sub is_bandpass_calibrator {
    my ($dbh,$source,$frequency)=@_;

    my $wavelength=&which_wavelength($frequency);
    my $isbandpass=-1;
    my $chk_query="SELECT bandpass_".$wavelength." FROM calibrators WHERE name='".$source."'";
    if ($output_debug>1){
	print "executing query: $chk_query\n";
    }
    my $chk=$dbh->prepare($chk_query);
    $chk->execute;
    while (my @chk_response=$chk->fetchrow){
	$isbandpass=$chk_response[0];
    }
    $chk->finish;

    return $isbandpass;
}

sub source_present {
    my ($source,@source_list)=@_;
    
    my $is_present=0;
    for (my $i=0;$i<=$#source_list;$i++){
	if ($source_list[$i] eq $source){
	    $is_present=1;
	    last;
	}
    }

    return $is_present;
}

sub add_source {
    my ($new_source,@source_list)=@_;

    my $is_new=1;
    for (my $i=0;$i<=$#source_list;$i++){
	if ($source_list[$i] eq $new_source){
	    $is_new=0;
	    last;
	}
    }

    if ($is_new==1){
	push @source_list,$new_source;
    }

    return @source_list;

}

sub array_to_text {
    # turns an array into a postgresql compliant input
    my (@array)=@_;

    my $string="{";
    # check for letters, but we don't look for 'E' which may be used
    # in exponential numbers
    my $needsquotes=0;
    for (my $i=0;$i<=$#array;$i++){
	if ($array[$i]=~/[a-zA-DF-Z]/){
	    $needsquotes=1;
	}
    }
    for (my $i=0;$i<=$#array;$i++){
	if ($i>0){
	    $string.=",";
	}
	if ($needsquotes==1){
	    $string.="'";
	}
	if (!$array[$i]){
	    if ($needsquotes==1){
		$string.="NA";
	    } else {
		$string.="999";
	    }
	} else {
	    $string.=$array[$i];
	}
	if ($needsquotes==1){
	    $string.="'";
	}
    }
    $string.="}";
    return $string;
}

sub get_id_info {
    my ($dbh,$get_id)=@_;

    my $get_query="SELECT * FROM observation WHERE id=".$get_id;
    my $get=$dbh->prepare($get_query);
    $get->execute;
    my @get_response;
    while (@get_response=$get->fetchrow){
	last;
    }
    $get->finish;

    return @get_response;
}

sub time_to_unixtime {
    my ($ptime)=@_;
    
    $ptime=~/^(\d{4})\-(\d{2})\-(\d{2})\s(\d{2})\:(\d{2})\:(\d{2})/;
    my ($year,$month,$day,$hour,$minute,$second)=($1,$2,$3,$4,$5,$6);
    $year-=1900;
    $month-=1;

    my $utime=mktime($second,$minute,$hour,$day,$month,$year,0,0);

    return $utime;
}

sub which_wavelength {
    my ($frequency)=@_;

    my @highlow_21cm=(0,2000);
    my @highlow_13cm=(2000,4000);
    my @highlow_6cm=(4000,6500);
    my @highlow_3cm=(6500,10000);
    my @highlow_12mm=(15000,25000);
    my @highlow_7mm=(25000,55000);
    my @highlow_3mm=(80000,110000);

    my $wavelength;
    if (($frequency>=$highlow_21cm[0])&&
	($frequency<$highlow_21cm[1])){
	$wavelength="21cm";
    } elsif (($frequency>=$highlow_13cm[0])&&
	     ($frequency<$highlow_13cm[1])){
	$wavelength="13cm";
    } elsif (($frequency>=$highlow_6cm[0])&&
	     ($frequency<$highlow_6cm[1])){
	$wavelength="6cm";
    } elsif (($frequency>=$highlow_3cm[0])&&
	     ($frequency<$highlow_3cm[1])){
	$wavelength="3cm";
    } elsif (($frequency>=$highlow_12mm[0])&&
	     ($frequency<$highlow_12mm[1])){
	$wavelength="12mm";
    } elsif (($frequency>=$highlow_7mm[0])&&
	     ($frequency<$highlow_7mm[1])){
	$wavelength="7mm";
    } elsif (($frequency>=$highlow_3mm[0])&&
	     ($frequency<$highlow_3mm[1])){
	$wavelength="3mm";
    }

    return $wavelength;
}

sub miriad_split_data {
    my ($atlod_data,$keep_data)=@_;

    my $uvsplit_command="uvsplit vis=".$atlod_data." options=clobber";
    if ($output_debug>1){
	print "MESSAGE: executing ".$uvsplit_command."\n";
    }
    my @uvsplit_output=&run_command($uvsplit_command);

    &print_output_lines("MIRIAD",@uvsplit_output);

    if ($keep_data!=1){
	# now remove the raw data, we won't be needing it anymore
	my $remove_command="rm -rf ".$atlod_data;
	if ($output_debug>1){
	    print "MESSAGE: executing ".$remove_command."\n";
	}
	system $remove_command;
    }

}

sub print_output_lines {
    my ($line_prefix,@lines)=@_;
    
    if ($output_debug>1){
	for (my $i=0;$i<=$#lines;$i++){
	    print $line_prefix.": ".$lines[$i]."\n";
	}
    }
}

sub is_planet {
    my ($source)=@_;

    $source=~tr/A-Z/a-z/; # change to lowercase

    my $planet=0;
    my @planets=qw/mercury venus earth mars jupiter saturn uranus neptune/;

    for (my $i=0;$i<=$#planets;$i++){
	if ($source eq $planets[$i]){
	    $planet=1;
	    last;
	}
    }

    return $planet;
}

sub miriad_calibrate {
    my ($source,$frequency,$refant,$do_bandpass,$do_gains,$fluxcal)=@_;

    # if the source is a planet, we don't do any calibration
    if (&is_planet($source)==1){
	if ($output_debug>1){
	    print "MESSAGE: calibration aborted; source is a planet\n";
	}
	return;
    }

    if ($do_bandpass==1){
	# we solve for the bandpass using mfcal
	my $mfcal_command="mfcal vis=".$source.".".$frequency.
	    " refant=".$refant." interval=0.1";
	if ($output_debug>1){
	    print "MESSAGE: executing ".$mfcal_command."\n";
	}
	my @mfcal_output=&run_command($mfcal_command);
	&print_output_lines("MIRIAD",@mfcal_output);
    }
    if ($do_gains==1){
	# we solve for the gains using gpcal
	my $calibrated=0;
	my $use_qusolve=1;
	while ($calibrated==0){
	    $calibrated=1;
	    my $gpcal_command="gpcal vis=".$source.".".$frequency.
		" refant=".$refant." interval=0.1 options=xyvary";
	    if (($fluxcal==0)&&($use_qusolve==1)){
		$gpcal_command.=",qusolve";
	    }
	    if ($output_debug>1){
		print "MESSAGE: executing ".$gpcal_command."\n";
	    }
	    my @gpcal_output=&run_command($gpcal_command);
	    &print_output_lines("MIRIAD",@gpcal_output);
	    # check that gpcal didn't fail
	    for (my $i=0;$i<=$#gpcal_output;$i++){
		if ($gpcal_output[$i]=~/^\#\#\#\sFatal\sError:.*degenerate$/){
		    $calibrated=0;
		    $use_qusolve=0;
		}
	    }
	}
    }
}

sub determine_reference_antenna {
    my ($dbh,@rpfits_files)=@_;

    my $ref_query="SELECT onsource FROM observation WHERE (";
    for (my $i=0;$i<=$#rpfits_files;$i++){
	if ($i!=0){
	    $ref_query.=" OR ";
	}
	$ref_query.="contains_text(rpfits_file,'".$rpfits_files[$i]."')=true";
    }
    $ref_query.=")";
    if ($output_debug>1){
	print "MESSAGE: executing query: $ref_query\n";
    }
    my $ref=$dbh->prepare($ref_query);
    $ref->execute;
    my @total_onsource_cycles;
    my $initialised_count=0;
    while (my @ref_response=$ref->fetchrow){
	my @onsource_cycles=&postgresql_array_get($ref_response[0]);
	if ($initialised_count==0){
	    for (my $i=0;$i<=$#onsource_cycles;$i++){
		$total_onsource_cycles[$i]=$onsource_cycles[$i];
	    }
	    $initialised_count=1;
	} else {
	    for (my $i=0;$i<=$#onsource_cycles;$i++){
		$total_onsource_cycles[$i]+=$onsource_cycles[$i];
	    }
	}	    
    }
    $ref->finish;
    my ($minants,$maxants)=&minmax(@total_onsource_cycles);
    my @order_preference=(2,3,1,4,5); # never allow CA06 to be reference
    my $reference_antenna=-1;
    while($reference_antenna==-1){
	for (my $i=0;$i<=$#order_preference;$i++){
	    if ($total_onsource_cycles[$order_preference[$i]-1]==$maxants){
		# this antenna will work as the reference
		$reference_antenna=$order_preference[$i];
		last;
	    }
	}
	$maxants--;
    }
    
    return $reference_antenna;
}

sub miriad_copy_calibration {
    my ($calibrated_source,$uncalibrated_source,$frequency)=@_;

    my $gpcopy_command="gpcopy vis=".$calibrated_source.".".
	$frequency." out=".$uncalibrated_source.".".$frequency;
    if ($output_debug>1){
	print "MESSAGE: executing ".$gpcopy_command."\n";
    }
    my @gpcopy_output=&run_command($gpcopy_command);
    &print_output_lines("MIRIAD",@gpcopy_output);
}

sub autoflag {
    my ($source,$frequency)=@_;

    my $miriad_mirflag_command="mirflag vis=".$source.".".$frequency.
	" options=amp,medsed,short";
    for (my $j=0;$j<2;$j++){
	my $mirflag_command=$miriad_mirflag_command;
	if ($j==0){
	    $mirflag_command.=" stokes=xx,xx,yy";
	} elsif ($j==1){
	    $mirflag_command.=" stokes=yy,xx,yy";
	}
	if ($output_debug>1){
	    print "MESSAGE: executing $mirflag_command\n";
	}
	my @mirflag_output=&run_command($mirflag_command);
	&print_output_lines("MIRIAD",@mirflag_output);
    }
}

sub event_flag {
    my ($dbh,$uvdata,@rpfits_files)=@_;

    my @event_times;

    # find the earliest and latest times available in the RPFITS files
    my $tim_query="SELECT time FROM observation WHERE (";
    for (my $i=0;$i<=$#rpfits_files;$i++){
	if ($i!=0){
	    $tim_query.=" OR ";
	}
	$tim_query.="contains_text(rpfits_file,'".$rpfits_files[$i]."')=true";
    }
    $tim_query.=") ORDER BY time LIMIT 1";
    if ($output_debug>1){
	print "MESSAGE: executing query: $tim_query\n";
    }
    my $tim=$dbh->prepare($tim_query);
    $tim->execute;
    my $earliest_time;
    while (my @tim_response=$tim->fetchrow){
	$earliest_time=$tim_response[0];
	last;
    }
    $tim->finish;
    if ($output_debug>0){
	print "MESSAGE: dataset starts at time ".$earliest_time."\n";
    }

    $tim_query="SELECT time+totaltime(maxcycles(onsource),maxcycles(offsource),".
	"cycle_time) FROM observation WHERE (";
    for (my $i=0;$i<=$#rpfits_files;$i++){
	if ($i!=0){
	    $tim_query.=" OR ";
	}
	$tim_query.="contains_text(rpfits_file,'".$rpfits_files[$i]."')=true";
    }
    $tim_query.=") ORDER BY time DESC LIMIT 1";
    if ($output_debug>1){
	print "MESSAGE: executing query: $tim_query\n";
    }
    $tim=$dbh->prepare($tim_query);
    $tim->execute;
    my $latest_time;
    while (my @tim_response=$tim->fetchrow){
	$latest_time=$tim_response[0];
	last;
    }
    $tim->finish;
    if ($output_debug>0){
	print "MESSAGE: dataset ends at time ".$latest_time."\n";
    }

    push @event_times,$earliest_time;
    # get the events that occurred during the files
    my $evs_query="SELECT * FROM events WHERE event_time>'".$earliest_time."' AND ".
	"event_time<'".$latest_time."' AND (event_type='dcal' OR ".
	"event_type='acal') ORDER BY event_time";
    if ($output_debug>1){
	print "MESSAGE: executing query: $evs_query\n";
    }
    my $evs=$dbh->prepare($evs_query);
    $evs->execute;
    while(my @evs_response=$evs->fetchrow){
	my $event_time=$evs_response[1];
	push @event_times,$event_time;
    }
    $evs->finish;
    push @event_times,$latest_time;

    # we now have a list of when each significant event occurred. we go through
    # this and determine the time difference between each event
    my $maxtimediff=0;
    my ($early_good_time,$late_good_time);
    for (my $i=0;$i<$#event_times;$i++){
	my $this_timediff=&time_difference($dbh,$event_times[$i],
					   $event_times[$i+1]);
	if (&time_bigger($dbh,$maxtimediff,$this_timediff)==1){
	    $maxtimediff=$this_timediff;
	    $early_good_time=$event_times[$i];
	    $late_good_time=$event_times[$i+1];
	}
    }
    # add a buffer to the earliest good time for safety
    my $cycle_time=&get_cycle_time($dbh,@rpfits_files);
    my $buffer_time=4*$cycle_time; # in seconds
    my $early_good_time=&time_add($dbh,$early_good_time,$buffer_time);

    if ($output_debug>0){
	print "MESSAGE: biggest block of good data goes for ".$maxtimediff." between ".
	    $early_good_time." and ".$late_good_time."\n";
    }

    # now flag around this good time block
    if ($early_good_time ne $earliest_time){
	my $uvflag_command="uvflag vis=".$uvdata." flagval=flag \"select=time(".
	    &postgresql_to_miriad_time($earliest_time).",".
	    &postgresql_to_miriad_time($early_good_time).")\"";
	if ($output_debug>1){
	    print "MESSAGE: executing ".$uvflag_command."\n";
	}
	my @uvflag_response=&run_command($uvflag_command);
	&print_output_lines("MIRIAD",@uvflag_response);
    }

    if ($late_good_time ne $latest_time){
	my $uvflag_command="uvflag vis=".$uvdata." flagval=flag \"select=time(".
	    &postgresql_to_miriad_time($late_good_time).",".
	    &postgresql_to_miriad_time($latest_time).")\"";
	if ($output_debug>1){
	    print "MESSAGE: executing ".$uvflag_command."\n";
	}
	my @uvflag_response=&run_command($uvflag_command);
	&print_output_lines("MIRIAD",@uvflag_response);
    }	

}

sub time_difference {
    my ($dbh,$early_time,$late_time)=@_;

    my $tdf_query="SELECT TIMESTAMP '".$late_time."' - TIMESTAMP '".$early_time."'";
    my $tdf=$dbh->prepare($tdf_query);
    $tdf->execute;
    my $timediff;
    while (my @tdf_response=$tdf->fetchrow){
	$timediff=$tdf_response[0];
	last;
    }
    $tdf->finish;

    return $timediff;
}

sub time_bigger {
    my ($dbh,$small_time,$big_time)=@_;

    my $tbg_query="SELECT INTERVAL '".$big_time."' > INTERVAL '".$small_time."'";
    my $tbg=$dbh->prepare($tbg_query);
    $tbg->execute;
    my $result;
    while (my @tbg_response=$tbg->fetchrow){
	$result=$tbg_response[0];
	last;
    }
    $tbg->finish;

    return $result;
}

sub postgresql_to_miriad_time {
    my ($postgresql_time)=@_;
    
    $postgresql_time=~/^([^\-]*)\-([^\-]*)\-([^\-]*)\s([^\:]*)\:([^\:]*)\:(.*)$/;
    my ($year,$month,$day,$hour,$minute,$second)=($1,$2,$3,$4,$5,$6);
    
    my $miriad_time=sprintf "%4d%3s%02d:%02d:%02d:%02d",$year,
    &month_number_to_name($month),$day,$hour,$minute,$second;
    
    return $miriad_time;
}

sub month_number_to_name {
    my ($month_number)=@_;
    
    my @month_names=("JAN","FEB","MAR","APR","MAY","JUN",
		     "JUL","AUG","SEP","OCT","NOV","DEC");
    $month_number--;
    return $month_names[$month_number];
    
}

sub time_add {
    my ($dbh,$reference_time,$interval)=@_;
    
    my $gtq_query="SELECT TIMESTAMP '".$reference_time."' + INTERVAL '".$interval." second'";
    if ($output_debug>1){
	print "MESSAGE: executing query: $gtq_query\n";
    }
    my $gtq=$dbh->prepare($gtq_query);
    $gtq->execute;
    
    my $final_time;
    while(my @times=$gtq->fetchrow){
	$final_time=$times[0];
    }
    $gtq->finish;

    return $final_time;
}

sub get_cycle_time {
    my ($dbh,@rpfits_files)=@_;
    
    my $gct_query="SELECT cycle_time FROM observation WHERE rpfits_file='".
	$rpfits_files[0]."' LIMIT 1";
    my $gct=$dbh->prepare($gct_query);
    $gct->execute;
    my $cycletime;
    while (my @gct_response=$gct->fetchrow){
	$cycletime=$gct_response[0];
	last;
    }
    $gct->finish;

    return $cycletime;
}

sub miriad_self_calibrate {
    my ($image_model,$source,$frequency,$maximum_flux,$refant)=@_;

    my $gpscal_command="gpscal vis=".$source.".".$frequency.
	" model=".$image_model." clip=".(0.01*$maximum_flux).
	" interval=1 options=amplitude,mfs,noscale,xyvary ".
	" refant=".$refant;
    if ($output_debug>1){
	print "MESSAGE: executing $gpscal_command\n";
    }
    my @gpscal_output=&run_command($gpscal_command);
    &print_output_lines("MIRIAD",@gpscal_output);
}

sub miriad_fluxscale {
    my ($fluxcal,$phasecal,$frequency)=@_;

    # is the flux calibrator a planet?
    my $fluxplanet=&is_planet($fluxcal);

    # set the gain scaling to be the same
    my $gpboot_command="gpboot vis=";
    if ($fluxplanet==0){
	# flux calibrator should have the correct amplitude
	$gpboot_command.=$phasecal.".".$frequency." cal=".$fluxcal.".".$frequency;
    } else {
	# flux calibrator will not have been gain calibrated
	$gpboot_command.=$fluxcal.".".$frequency." cal=".$phasecal.".".$frequency;
    }
    if ($output_debug>1){
	print "MESSAGE: executing ".$gpboot_command."\n";
    }
    my @gpboot_output=&run_command($gpboot_command);
    &print_output_lines("MIRIAD",@gpboot_output);

    # now correct the spectral slope
    my $mfboot_command="mfboot vis=".$fluxcal.".".$frequency.",".
	$phasecal.".".$frequency." \"select=source(".$fluxcal.")\"";
    if ($output_debug>1){
	print "MESSAGE: executing ".$mfboot_command."\n";
    }
    my @mfboot_output=&run_command($mfboot_command);
    &print_output_lines("MIRIAD",@mfboot_output);
    
}

sub miriad_make_image {
    my ($imaging_source,$xpixels,$ypixels,$frequency,$robust_value)=@_;

    my $imagetripled=0;

    my $dataprefix=$imaging_source.".".$frequency;
    my $map=$dataprefix.".imap";
    my $beam=$dataprefix.".ibeam";
    # check if the datasets already exist, and delete them if they do
    if (-e $map){
	my $remove_command="rm -rf ".$map;
	if ($output_debug>1){
	    print "MESSAGE: executing ".$remove_command."\n";
	}
	system $remove_command;
    }
    if (-e $beam){
	my $remove_command="rm -rf ".$beam;
	if ($output_debug>1){
	    print "MESSAGE: executing ".$remove_command."\n";
	}
	system $remove_command;
    }

    my $invert_command="invert vis=".$dataprefix.
	" map=".$map." beam=".$beam." stokes=i robust=".$robust_value;
#    my ($xpixels,$ypixels)=(-1,-1);
#    if (($xbeamsize!=-1)&&($ybeamsize!=-1)){
#	# we have information about the size of the synthesised beam, so
#	# we can make a nicer image with square pixels and use the sdb option
#	# as well
#	if ($output_debug>1){
#	    print "MESSAGE: beamsize is ".$xbeamsize." x ".$ybeamsize.
#		" arcsecs\n";
#	}
#	my $xcell=$xbeamsize/3.0;
#	my $ycell=$ybeamsize/3.0;
#	$invert_command.=" cell=".$xcell.",".$ycell;
#	# compute the size of the primary beam
#	my $pbeamsize=&primary_beamsize($frequency); # comes back as radians
#	$pbeamsize*=180.0/3.141592654; # now in degrees
#	$pbeamsize*=3600.0; # now in arcseconds, same as the cell size
#	$xpixels=floor($pbeamsize/$xcell);
#	$ypixels=floor($pbeamsize/$ycell);
#	# since we're going to make an MFS image with a SDB, so we can use
#	# mfclean, we'll make the image 3 times bigger than the primary beam
#	$xpixels*=3;
#	$ypixels*=3;
    if (($xpixels!=-1)&&($ypixels!=-1)){
	$xpixels*=3;
	$ypixels*=3;
	$invert_command.=" imsize=".$xpixels.",".$ypixels;
	$invert_command.=" options=systemp,mfs,sdb";
	$imagetripled=1;
    } else {
	# this is probably the first time we've been asked to image, so
	# we let MIRIAD do what it thinks is best
	$invert_command.=" options=double,systemp,mfs";
    }

    if ($output_debug>1){
	print "MESSAGE: executing ".$invert_command."\n";
    }
    my @invert_output=&run_command($invert_command);
    &print_output_lines("MIRIAD",@invert_output);

    # find the theoretical noise of the image
    my $theoretical_rms=-1;
    for (my $i=0;$i<=$#invert_output;$i++){
	if ($invert_output[$i]=~/^Theoretical\srms\snoise\:\s(.*)$/){
	    $theoretical_rms=$1;
	}
    }

    # return with the name of the image, beam and theoretical noise, along
    # with the number of pixels in each image direction
    return ($map,$beam,$theoretical_rms,$xpixels,$ypixels,$imagetripled);
}

sub miriad_apply_calibration {
    my ($source,$frequency)=@_;

    # use uvaver to apply the calibration tables to the visibility data
    my $uvaver_command="uvaver vis=".$source.".".$frequency.
	" out=".$source.".".$frequency.".uvaver";
    if ($output_debug>1){
	print "MESSAGE: executing ".$uvaver_command."\n";
    }
    my @uvaver_output=&run_command($uvaver_command);
    &print_output_lines("MIRIAD",@uvaver_output);

    # write the applied dataset over the unapplied one
    my $remove_command="rm -rf ".$source.".".$frequency;
    if ($output_debug>1){
	print "MESSAGE: executing ".$remove_command."\n";
    }
    system $remove_command;
    my $rename_command="mv ".$source.".".$frequency.".uvaver ".
	$source.".".$frequency;
    if ($output_debug>1){
	print "MESSAGE: executing ".$rename_command."\n";
    }
    system $rename_command;

}

sub miriad_deconvolve_image {
    my ($source,$frequency,$map,$beam,$rms,$xpixels,$ypixels)=@_;

    my $dataprefix=$source.".".$frequency;
    my $model=$dataprefix.".iclean";
    if (-e $model){
	my $remove_command="rm -rf ".$model;
	if ($output_debug>0){
	    print "MESSAGE: executing ".$remove_command."\n";
	}
	system $remove_command;
    }
    my $clean_command;
    if (($xpixels==-1)&&($ypixels==-1)){
	# no idea about the image size, and invert was used to make a
	# normal dirty map, so the beam has only one plane, and we have
	# to use regular clean
	$clean_command="clean map=".$map." beam=".$beam." out=".$model.
	    " cutoff=".(5*$rms)." speed=-0.5 niters=3000 options=positive";
    } else {
	# invert was used to make a spectral dirty beam, so we use
	# mfclean
	$clean_command="mfclean map=".$map." beam=".$beam." out=".$model.
	    " cutoff=".(5*$rms)." speed=-0.5 niters=3000 ";
	my $halfx=floor($xpixels/6.0);
	my $halfy=floor($ypixels/6.0);
	$clean_command.="\"region=relcenter,boxes(-".$halfx.",-".$halfy.",".
	    $halfx.",".$halfy.")\"";
    }
    if ($output_debug>0){
	print "MESSAGE: executing ".$clean_command."\n";
    }
    my @clean_output=&run_command($clean_command);
    &print_output_lines("MIRIAD",@clean_output);

    # return the name of the model
    return $model;
}

sub miriad_restore_image {
    my ($source,$frequency,$map,$beam,$model)=@_;

    my $dataprefix=$source.".".$frequency;
    my $image=$dataprefix.".irestor";
    if (-e $image){
	my $remove_command="rm -rf ".$image;
	if ($output_debug>0){
	    print "MESSAGE: executing ".$remove_command."\n";
	}
	system $remove_command;
    }
    my $restor_command="restor model=".$model." beam=".$beam.
	" map=".$map." out=".$image;
    if ($output_debug>0){
	print "MESSAGE: executing ".$restor_command."\n";
    }
    my @restor_output=&run_command($restor_command);
    &print_output_lines("MIRIAD",@restor_output);

    # look for the size of the synthesised beam
    my ($semimajor,$semiminor,$positionangle);
    for (my $i=0;$i<=$#restor_output;$i++){
	if ($restor_output[$i]=~/^Using\sgaussian\sbeam\sfwhm\sof\s+(.*)\sby\s+(.*)\sarcsec\.$/){
	    $semimajor=$1;
	    $semiminor=$2;
	} elsif ($restor_output[$i]=~/^Position\sangle\:\s+(.*)\sdegrees\.$/){
	    $positionangle=$1;
	}
    }

    # calculate the synthesised beam size in the RA and Dec directions
    # remembering that position is measured east from north
    my $maj_dir=0; # by default major direction is along Declination
    my $min_dir=1; # minor direction is along RA
    if ($positionangle<0){
	$positionangle+=360.0;
    }
    while($positionangle>90){
	if ($maj_dir==0){
	    $maj_dir=1;
	    $min_dir=0;
	} else {
	    $maj_dir=0;
	    $min_dir=1;
	}
	$positionangle-=90.0;
    }
    my ($ra_axis,$dec_axis);
    my $deg_to_rad=3.141592654/180.0;
    if ($maj_dir==0){
	$dec_axis=$semimajor*cos($positionangle*$deg_to_rad);
	$ra_axis=$semiminor*cos($positionangle*$deg_to_rad);
    } else {
	$ra_axis=$semimajor*cos($positionangle*$deg_to_rad);
	$dec_axis=$semiminor*cos($positionangle*$deg_to_rad);
    }

    # return the size of the synthesised beam and the name of the image
    return ($image,$ra_axis,$dec_axis);
}

sub miriad_image_size {
    my ($image)=@_;

    my ($xpixels,$ypixels,$max_flux);

    my $imlist_command="imlist in=".$image." options=statistics";
    my @imlist_output=&run_command($imlist_command);
    &print_output_lines("MIRIAD",@imlist_output);

    my $reached_planes=0;
    $max_flux=0;
    for (my $i=0;$i<=$#imlist_output;$i++){
	my $line=$imlist_output[$i];
	if ($line=~/Trc\=\((.*)\,(.*)\,(.*)\)/){
	    $xpixels=$1;
	    $ypixels=$2;
	} elsif ($line=~/plane\s+Frequency/){
	    $reached_planes=1;
	} elsif ($reached_planes==1){
	    $line=~s/^\s*//;
	    my @els=split(/\s+/,$line);
	    if ($els[3]>$max_flux){
		$max_flux=$els[3];
	    }
	}
    }

    return ($xpixels,$ypixels,$max_flux);
}

sub primary_beamsize {
    my ($frequency)=@_;

    my $dishdiameter=22; # metres, the size of each ATCA antenna
    my $c=299792458; # m/s, the speed of light
    $frequency*=1e6; # Hz, the frequency of the observation

    my $wavelength=$c/$frequency; # m, the wavelength of the observation

    my $beamsize=$wavelength/$dishdiameter; # radians
}

sub obtain_rpfits {
    my ($rpfits_file,$work_directory)=@_;

    # try to find the RPFITS file
    my @search_paths=("/DATA/ARCHIVE_1/ATCA/archive");
    my $found_file="";
    $rpfits_file=~s/^.*\/(.*)$/$1/;

    for (my $i=0;$i<=$#search_paths;$i++){
	$search_paths[$i]=~s/\/$//;
	my $search_file=$search_paths[$i]."/".$rpfits_file;
	if (-e $search_file){
	    $found_file=$search_file;
	    last;
	}
    }
    if ($found_file eq ""){
	# couldn't find this file!
	if ($output_debug>1){
	    print "ERROR: unable to find RPFITS file $rpfits_file from local repository!\n";
	}

	# try to get this file from the ATOA
	$rpfits_file=~/(\d{4})\-(\d{2})\-(\d{2})\_(.*)/;
	my $rpfits_year=$1;
	my $rpfits_month=$2;
	my @rpfits_archive_search=("/data/ATOA_1/ATCA",
				   "/data/ATOA_2/ATCA");
	for (my $i=0;$i<=$#rpfits_archive_search;$i++){
	    my $rpfits_archive_path=$rpfits_archive_search[$i];
	    $rpfits_archive_path.="/".$rpfits_year."-".$rpfits_month;

	    my $copy_command="rsync -av ste616\@horus:".$rpfits_archive_path."/";
	    $copy_command.=$rpfits_file;
	    $copy_command.=" ".$work_directory."/";
	    if ($output_debug>0){
		print "executing: $copy_command\n";
	    }
	    system $copy_command;
	    
	    $found_file=$work_directory."/".$rpfits_file;
	    if (!-e $found_file){
		$found_file="";
	    } else {
		last;
	    }
	}
    }

    if ($found_file eq ""){
	# still couldn't find the file
	if ($output_debug>0){
	    print "ERROR: unable to find RPFITS file $rpfits_file from any repository!\n";
	}	
	return "";
    }

    # return the location of the file we just extracted
    return $found_file;
    
}
